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ABSTRACT 

Since the late 1980’s, research in recursive formulations of multibody dynamics has flourished. 
Historically, much of this research can be traced to applications of low dimensionality in mecha- 
nism and vehicle dynamics. Indeed, there is littfe doubt that recursive order N methods are the 
method of choice for this class of systems. This approach has the advantage that a minimal num- 
ber of coordinates are utilized, parallelism can be induced for certain system topologies, and the 
method is of order N computational cost for systems of N rigid bodies. 

Despite the fact that many authors have dismissed redundant coordinate formulations as being of 
order N 3 , and hence less attractive than recursive formulations, we present recent research that 
demonstrates that at least three distinct classes of redundant, nonrecursive multibody formulations 
consistently achieve order N computational cost for systems of rigid and/or flexible bodies. These 
formulations are the 

(i) preconditioned range space formulation, 

(ii) penalty methods, and 

(iii) augmented Lagrangian methods 

for nonlinear multibody dynamics. The first method can be traced to its foundation in equality 
constrained quadratic optimization, while the last two methods have been studied extensively in 
the context of coercive variational boundary value problems in computational mechanics. Until 
recently, however, they have not be investigated in the contrext of multibody simulation, and 
present theoretical questions unique to nonlinear dynamics. All of these nonrecursive methods 
have additional advantages with respect to recursive order N methods: (1) The formalisms retain 
the highly desirable order N computational cost, (2) the techniques are amenable to concurrent 
simulation strategies, (3) the approaches do not depend upon system topology to induce concur- 
rency, (4) and the methods can be derived to balance the computational load automatically on con- 
current multiprocessors. In addition to the presentation of the fundamental formulations, this 
paper presents new theoretical results regarding the rate of convergence of order N constraint sta- 
bilization schemes associated with the newly introduced class of methods. 
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Introduction 


It is well known that the nonlinear equations governing the motion of complex multibody systems 
give rise to differential-algebraic equations that present unique and interesting problems in their 
solution [Gear, Petzold], The mathematical study of differential-algebraic equationsis justifiably a 
field of research unto itself. However, in multibody applications, it appears that most researchers 
prefer to eliminate the troublesome multipliers, and deal exclusively with the problems associated 
with the solution of sets of ordinary differential equations. The point in the analysis at which the 
multipliers are eliminated has become one criteria for distinguishing among the various formula- 
tions. During the last 1970’s and early 1980’s, numerous publications appeared in which the con- 
straint multipliers are eliminated at computation time by numerically constructing a basis for the 
instantaneous nullspace of the constraint Jacobian. Some representative work from this class 
includes [Singh, Wampler, Wehage, Amirouche,...], and are referred to as the nullspace methods. 
These algorithms have been so named due to their similarity to the nullspace methods in qua- 
dratic, equality constrained optimizations [Gill]. Despite their elegance, numerical calculation of 
the nullspace basis and its use to eliminate the multipliers leads to a dense system coefficient 
matrix. This matrix is of order (N-M) X (N-M) where N is the number of redundant coordinates 
and M is the number of constraints. Consequently, the nullspace methods are of O(N-M) 3 at each 
time step. 

As opposed to techniques that eliminate the multipliers numerically at computation time, a number 
of authors have derived elegant techniques for eliminating the constraints a priori; that is during the 
derivation of the equations. The works [Rodriguez, June 1987; Rodriguez, 1987] are 
representative of this class. These methods have a number of distinct and well-known advantages 
over the cubic order nullspace methods: 

(i) They employ a minimal coordinate set. 

(ii) They attain an O(N) computational cost for systems of rigid bodies. 

(iii) They can be employed in concurrent architectures for classes of problems. 

These methods have come to be known as the order N or recursive methods of formulating multi- 
body dynamics. 

Despite their numerous advantages, there are several aspects of the recursive order N methods that 
can cause difficulties for classes of problems. Foremost among these problems is that concurrency 
in the recursive methods is induced, at present, by assigning topologically independent branches 
of the structure to computationally independent processors. One need only consider the problem 
of modelling a tethered satellite system to realize that not all structures exhibit such structural par- 
allelism. Furthermore, most systems exhibit low levels of inherent structural parallelism appropri- 
ate for the concurrent implementation of recursive order N methods as in [Bae]. Typical space 
structures such as the space station or the CSI Evolutionary model at NASA Langley have only on 
the order of 6-12 “independent branches.” Gross underutilization of moderately parallel architec- 
tures, with 16 to 128 processors, can result when mapping the recursive algorithm onto such a 
computer. Finally, computational load balancing of the problem among processors in most imple- 
mentation is carried out by hand, which can be a tedious and time consuming task. 
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Nonrecursive Formulations 

This paper surveys the results of several recent approaches presented in [Kurdila 1,2,3] and 
described in detail in the dissertation by [Menon]. Collectively, the authors refer to these 
methods as nonrecursive formulations of multibody dynamics. In point of fact, the nonrecursive 
formulations described in this paper represent a synthesis of a family of formulations and con- 
straint stabilization procedures that are closely related. These methods can be employed simulta- 
neously independently, or in combination selected by an analyst to meet accuracy and 
computational speed requirements in a highly predictable manner. This class includes 

(j) The Prr./Ranpe Spare Method : As with the nullspace method, this 
approach can be traced historically to equality constrained quadratic optimization 
[Gill]. Analysis of the application of the approach to multibody simulation has been 
carried out recently in [Kurdila] and [Menon] using an iterative technique, and earlier 
noniterative versions appear in [Placek ] and [Wittenburg]. The primary advantage of 
the range space method is that is is the most numerically efficient of the three methods 

presented in this paper. 

(ii) Th r Penalty Method : Of course, the penalty method has been studied in detail 
in application to optimization [Luenberger] and coercive variational boundary value 
problems [Oden]. Still, it has only recently been studied in detail in the context of 
nonlinear multibody simulation [Bayo, Park, Kurdila]. An advantage of this approach 
is that explicit bounds on the constraint violation can be achieved in terms of the 
penalty parameter, even in the nonlinear case. [Kurdila] One disadvantage of the 
method that has been noted in the literature [Kurdila] , [Bayo] is that prohibitively 
large values of the penalty parameter are required to achieve the analytical guarantees 
of accuracy, but result in numerical ill-conditioning. 

(iii ) Thf Augmented T ^grangia n Method : As in the case of the penalty method, 
the augmented Lagrangian formulations have been studied extensively for use in 
the solution of variational boundary value problems [Glowinski]. Again, as in the 
case of the penalty method, this approach has only recently been studied within 
the field of nonlinear multibody dynamics [Bayo]. Empirical evidence in [Bayo 1,2] 
suggests that the method is superior to the penalty method from a computational 
viewpoint in that required accuracy can be achieved without large penalty parameters. 
The analytical study of Augmented Lagrangian formulations of nonlinear multibody 
dynamics remains an open and interesting field of research. 

To derive the equations employed in this paper, one must first consider the exact governing system 
of differential-algebraic equations 


M (q) q =/ (q> q,t)+C?X 

=0 


( 1 ) 
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and differentiated form of the constraints 


ao 

a q J 


q + ® t = 0 


( 2 ) 


C ( q , 0 


ao* 

.dq. 


(3) 


In these equations, q are the generalized coordinates selected for the problem at hand, and <t>(q) 
are the functional relations defining the holonomic constraints acting on the system.While the 
interested reader is referred to [Menon] for the details of their derivation, the equations govern- 
ing the simulation of multibody systems considered in this paper can be written 

X o = 0 

or 

(CAT l C T ) X 0 = -CAT'f-Cq-®' 


(4) 

(5) 


M (?) 4 n =/(«> t) -C t) - (<i> n + + to 2 <D) -cfx n 


1 


X n + l =X n + p (<i > „ + 2 ^0) < i> + to : Z <D) 


( 6 ) 

(7) 


The terms X are unknown Lagrange multipliers, X n+l are iterative corrections the the multipli- 
ers, and are iterative values of the constraint equations. The relationship of the above system 
to the aforementioned range space, penalty and Augmented Lagrangian formulations can be sum- 
marized as follows: 

(i) By employing equation (5) to initiate the integration procedure and letting 

l = 0 <8) 
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the above equations reduce to the range space equations, and provide the most computationally 
efficient version of the above set of equations when a judiciously selected preconditioned conju- 
gate gradient method is employed. 

(ii) By employing only equation (6), and by setting 

^ = X n+l S ° (9) 


the equations above reduce to the inertial penalty method introduced in [Bayo] and studied in fur- 
ther detail in [Kurdila]. This form of the governing equations provides some guarantees on 
accuracy, but is more expensive computationally. 

(iii) By using equation (4) for the initialization procedure, and choosing equations (6) 
and (7) iteratively, these equations reduce to the augmented Lagrangian approach introduce in 
[Bayo] and studied further in [Kurdila] and [Menon], While this version of the governing equa- 
tions are the most computationally expensive, they also provide the most control over the accu- 
racy attained in the simulation. 

At this point it should be noted that both the penalty and augmented Lagrangian methods can be 
viewed as stabilizaton procedures for the range space method. In the following discussion, the 
range space formulation will be referred to as the “baseline method”. In fact, it will be demon- 
strated later that this is a useful interpretation and leads to hybrid simulation methods that com* 
bine the computational efficiency of the range space method and the accuracty of the penalty or 
augmented Lagrangian methods. 

In the following sections it will be emphasized that the simulation of transient response of multi- 
body systems using various versions of the nonrecursive formulation above has many advantages: 

(i) The method achieves an order N computational cost while employing a 
redundant formulation. 

(ii) Consequently, the often complicated relative reference frame kinematics 
of the recursive methods can be avoided. 

(iii) The method is amenable to concurrent implementation on a wide variety of 
forthcoming parallel architectures. 

(iv) The technique is essentially automatically, computationally load balancing 

(v) The approach does not depend upon system topology to induce parallelism, and 
does not result in gross underutilization of available processors. 
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Order N Computational Cost of the Baseline Algorithm 

The order N computations] cost of the baseline algorithm when the governing equations are 
selected to be simplified versions of equations (5) and (6) 

( CAT 1 C T ) l = -CAT l f-Cq-& t 

. ( 10 ) 

M (q)q=f(q, q, t) ~(f\ 

, otherwise denoted as the preconditioned range space formulation in [Kurdila] and [Menon], is 
well documented. While the details of the convergence properties of the algorithm exceed the 
scoprof this paper the reader is referred to [Menon] for a complete discussion In summary, the 
computational performance of the algorithm can be attributed to 

(i) the derivation of a rapidly convergent block Jacobi preconditioner based on 
the directed connectivity graph of the multibody system, 

(ii) the exploitation of the block-diagonal structure of the system coefficient 
matrices in the formulation, and 

(iii) the parallel implementation of the algorithm based upon subdomain 
decomposition techniques that are only weakly coupled to the system topology. 

Results extracted from [Kurdila] and [Menon] in figures (1) through (4) illustrate results typical of 
the preconditioned conjugate gradient / range space formulation. 

Accuracy and Order N Performance with Constraint Stabilization 

The motivation for deriving additional variants of the preconditioned conjugate gradient / range 
space formulation arises from the well-known need for stabilization in redundant formulations as 
well as the documented conditioning problems that can occur in the range space method of opti- 
mization [Gill]. Essentially, the accuracy of the simulation relies on the condition number of the 
constraint metric 

CM~ x C t ( 11 ) 

appearing in equation (5) and (10). As noted in [Gilj^ this condition number is bounded above by 


k(CAT 1 C 7 ') <k 2 (C)k(M) 


( 12 ) 


and may become large when the constraints become nearly redundant. However, inasmuch as 
great efforts have been made to ensure the order N computational cost of the baseline formulation, 
it has been a central goal in this work to derive stabilization methods that retain the order N corn- 
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Figure (1) : Order N Timing Results for Z- Truss 

Prototype Model 
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Figure (2) : Order N Timing Results for Space 
Station Prototype Model 



192 







Figure (3) : Order N Timing Results for CSI Evolutionary 
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putational cost. 

This goal has in fact been achieved and is depicted in figures (5) and (6). In figure (5), the timing 
results for a benchmark linear structural problem as a function of the number of degrees of free- 
dom is depicted. Two critical observations should be made upon inspection of this graph 

(i) All variants of the class of nonrecursive algorithms under investigation exhibit an 
order N compuational cost. 

(ii) The baseline preconditioned conjugate gradient / range space algorithm is the most 
efficient simulation. 

(iii) The computational cost of the nonrecursive methods with stabilization increases 
with the number of iterative corrections per time step. Hence, the penalty stabilization 
is more expensive thant the preconditioned conjugate gradient / ranges space method. 
But, the augmented Lagrangian formulation is more computationally expensive than 
the penalty method. In fact, the slope of the performance curve simply rotates 
counterclockwise as the number of fixed iterations per time step increases. 

In view of the measures of accuracy depicted in figures (7) and (8), it is not surprising that the iter- 
ative corrections of the augmented Lagrangian are more costly. It is precisely the addition of the 
stabilization methods that yields the desired improvement in accuracy. As shown in [Kurdila], for 
certain classes of multibody systems the constraint violation can be bounded by 


| 0 > 


<J>J 2 < 


min (cr 


2E(0) 

, (a) ,a 2 . 

m v tnin 


(P)) 


As noted earlier, a major criticism of achieving stringent tolerances using this bound is that large 
penalty factors must be used in the simulation which can lead to poor conditioning in practice. 
This result has been extended rigorously to the augmented Lagrangian stabilization methods in 
[Bustimante] and [Menon] in the form of the equation 

|$ + 2^to6 + co 2 o| . + j £ + 2^toO + (0 2 o| ( . 


and is depicted graphically in figure (7). Furthermore, empirical results suggest the stronger result 
that 


<j>(t ) 2 + co 2 <I>(/) 2 | j + * 


2k 

< (-) |< j >(/) 2 + co 2 0 ( 0 2 

v \r 1 
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Sequential CPU time comparison for truss model 
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Sequential Vs Parallel Comparison for truss model 
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Simply put, the “energy of constraint violation shrinks by a factor of e for each iteration” 
employed in the augmented Lagrangian formulation. 

Thus, figures (5) and (7), illustrate the fundamental tradeoff between execution time and accuracy. 
Another natural consequence of this tradeoff is the design of hybrid methods such as depicted in 
figure (8). In these methods, the preconditioned conjugate gradient / range space method is 
empolyed until a specific pre-selected tolerance is met, at which point the augmented Lagrangian 
method is employed until another more stringent tolerance is reached. In this fashion, the compu- 
tational efficiency of the range space method is retained as long as possible, at which point the 
iterative augmented Lagrangian method is used to guarantee the desired accuracy. 

Conclusions 

This paper has presented a class of nonrecursive formulations of multibody dynamics, all of 
which exhibit an order N computational cost. The methods are complementary to the existing 
class of recursive order N methods in that they seem most appropriate for large dimensional mod- 
els without inherent structural parallelism. Figures (9) and (10) graphically summarize the contri- 
bution of this paper. Figure (9) shows that as of 1988 only one order N method had been derived. 
Figure (10) depicts the state of affairs currently, and shows that at least three distinct approaches 
can achieve order N computational cost, as well as numerous combinations of these algorithms. 
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